#include <stdio.h>
#include <conio.h>
#include <math.h>
#include <memory.h>

#define DBG_PRINT(_X_)  printf _X_
//#define DBG_PRINT(_X_)
#define EPSILON     0.000001
//#define EPSILON     0.000000001
#define MATRIX_SIZE 4




void gaus(double a[MATRIX_SIZE][MATRIX_SIZE],double *b, double *res)
{
    double c[MATRIX_SIZE],x[MATRIX_SIZE],wer;

    int i,j,l;
    DBG_PRINT(("\nOutput isxodnik\n"));

    for(i=0;i<MATRIX_SIZE;i++)
    {
        for(j=0;j<MATRIX_SIZE;j++)
        {
            DBG_PRINT(("\t%2.2f ",a[i][j]));
        };
        DBG_PRINT(("\t| %2.2f",b[i]));
        DBG_PRINT(("\n"));
    };

    for(l=0;l<MATRIX_SIZE;l++)
    {
        for(i=l;i<MATRIX_SIZE;i++)
        {
            c[i]=a[i+1][l]/a[l][l];
            b[i+1]=b[i+1]-b[l]*c[i];
            for(j=l;j<MATRIX_SIZE;j++)
            {
                a[i+1][j]=a[i+1][j]-a[l][j]*c[i];
            };
        };
    };

    DBG_PRINT(("\nTreugolnik:\n"));

    for(i=0;i<MATRIX_SIZE;i++)
    {
        for(j=0;j<MATRIX_SIZE;j++)
        {
            DBG_PRINT((" %2.2f ",a[i][j]));
        };
        DBG_PRINT((" | %2.2f",b[i]));
        DBG_PRINT(("\n"));
    };

    DBG_PRINT(("\n\n"));

    for(i=0;i<MATRIX_SIZE;i++)    
    {
        x[i]=0;
    };

    for(i=(MATRIX_SIZE-1);i>=0;i--)
    {
        wer=0;
        if(i!=4) for(j=i+1;j<MATRIX_SIZE;j++)
        {
            wer=wer+a[i][j]*x[j];
        };

        x[i]=(b[i]-wer)/a[i][i];
    };

    DBG_PRINT(("Reshenie:\n"));
    for(i=0;i<MATRIX_SIZE;i++)
    {
        DBG_PRINT(("x[%d] = %2.2f\n",i,x[i]));
        res[i]=x[i];
    };
}

void calculate_variant(double a[MATRIX_SIZE][MATRIX_SIZE],double b[MATRIX_SIZE])
{
    //double f[]={0.11,2.00,0.11,0.13};
    double f[MATRIX_SIZE];
    int i,j;
    double temp_a[MATRIX_SIZE][MATRIX_SIZE];
    double temp_b[MATRIX_SIZE];
    double temp_f[MATRIX_SIZE],f0[MATRIX_SIZE],xx[MATRIX_SIZE],delta_xx[MATRIX_SIZE];
    bool isExit;

    memcpy(f,b,sizeof(f));


    for(i=0;i<MATRIX_SIZE;i++)
    {
        for(j=0;j<MATRIX_SIZE;j++)
        {
            printf(" %2.2f ",a[i][j]);
        };
        printf(" | %2.2f",b[i]);
        printf("\n");
    };

    //----------------------

    memcpy(temp_b,b,sizeof(temp_b));
    memcpy(temp_a,a,sizeof(temp_a));
    gaus(temp_a,temp_b,xx);


    int count = 0;
    do
    {
        DBG_PRINT(("\nf0:\n"));
        for(i=0;i<MATRIX_SIZE;i++)                
        {              
            f0[i]=0;
            for(j=0;j<MATRIX_SIZE;j++)
            {
                f0[i]=f0[i]+(a[i][j]*xx[j]);
            };
            DBG_PRINT(("%2.5f ",f0[i]));
        };

        DBG_PRINT(("\nb:\n"));

        for(i=0;i<MATRIX_SIZE;i++) {
            temp_b[i]=(f0[i]-f[i])*(-1);
            DBG_PRINT(("%2.5f ",temp_b[i]));
        }

        memcpy(temp_a,a,sizeof(temp_a));
        gaus(temp_a,temp_b,delta_xx);

        DBG_PRINT(("\n delta_xx: \n"));
        for(i=0;i<MATRIX_SIZE;i++) {
            DBG_PRINT((" %2.5f\n",delta_xx[i]));
        }

        isExit = true;
        for(i=0;i<MATRIX_SIZE;i++){
            DBG_PRINT(("%2.5f ",delta_xx[i]));
            if (delta_xx[i]>EPSILON)
                isExit = false;
        }
        if (isExit == true)
            break;


        DBG_PRINT(("\ndelta:\n"));
        for(i=0;i<MATRIX_SIZE;i++){
            xx[i]=xx[i]+delta_xx[i];
        }
        if (isExit == true)
            break;

        for(i=0;i<MATRIX_SIZE;i++)
            f[i]=f0[i];
        count++;
    }
    while(count<200000);

    printf("\nPROVERKA______!!!!!!!!!!!!: count[%d]\n",count);
    for(i=0;i<MATRIX_SIZE;i++)                
    {                      
        f0[i]=0;
        for(j=0;j<MATRIX_SIZE;j++)
        {
            f0[i]=f0[i]+(a[i][j]*xx[j]);
        };
        printf("%2.2f ",f0[i]);
    };

    getch();
}


void calculate_lab1()
{
    
    double b[]={0.11,2.00,0.11,0.13};
    double a[4][4]=
    {
        {0.17,0.75,-0.18,0.21},
        {0.75,0.13,0.11,1.00},
        {-0.33,0.11,3.01,-2.01},
        {0.11,1.12,1.11,-1.31}
    };
    
    /*
    double b[]={11.33,32,42};
    double a[MATRIX_SIZE][MATRIX_SIZE]=
    {
        {6,-1,-1},
        {-1,6,-1},
        {-1,-1,6}
    };
    */
    calculate_variant(a,b);
}